Symbolic computation of second-order normal forms for 
Hamiltonian systems relative to periodic flows 

M. Avendano-Camacho, J. A. Vallejo and Yu. Vorobjev 
' Departamento de Matematicas, Universidad de Sonora (Mexico) 

CN ■ Blvd L. Encinas y Rosales s/n Col. Centro 

Edificio 3K-1 CP 83000 Hermosillo (Son) 
Email (JAV, corresponding author): jvallejo@fc.uaslp.mx 



m ! March 26, 2013 



Oh! 



Abstract 

A Maxima package called pdynamics is described. It is aimed to study Poisson (and 
symplectic) systems and, particularly, the determination of the second-order normal 
form for perturbed Hamiltonians H t — Hq + eH\ +e 2 i?2, relative to the periodic flow of 
the unperturbed Hamiltonian Hq. The formalism presented here is global, it does not 
require recursive computations and allows an efficient symbolic implementation. 
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1 Introduction 



In this paper we discuss some computational aspects of the normal form theory for Hamil- 
tonian systems on general phase spaces, that is, Poisson manifolds. According to Deprit 
[9], a perturbed vector field 

2 k 

A = A + eA 1 + e -A 2 + --- + e -A k + 0(e k+1 ) 

on a manifold M, is said to be in normal form of order k relative to Aq if p4o>-<4i] = 
for i € {1, . . . , k}. In the context of perturbation theory, the normalization problem is 
formulated as follows: to find a (formal or smooth) transformation which brings a perturbed 
dynamical system to a normal form up to a given order. 

The construction of a normalization transformation, in the framework of the Lie transform 
method [El [T2| [T3l 115] . is related to the solvability of a set of linear non homogeneous 
equations, called the homological equations. If the homological equations admit global 
solutions, defined on the whole M, we speak of a global normalization, which essentially 
depends on the properties of the unperturbed dynamics. 

Here we are interested in the global normalization of a perturbed Hamiltonian dynamics 
relative to periodic Hamiltonian flows. In this case, a result due to Cushman [6], states that 
if A is Hamiltonian, and the flow of the unperturbed vector field Aq is periodic, then the 
true dynamics admits a global Deprit normalization to arbitrary order. The corresponding 
normal forms can be determined by a recursive procedure (the so-called Deprit diagram) 
involving the resolution of the homological equations at each step. 

In this paper, we extend Cushman's result to the Poisson case and derive an alterna- 
tive coordinate-free representation for the second-order normal form, involving only three 
intrinsic operations: the averaging operators associated to the S 1 — action, and the Pois- 
son bracket. We give a direct derivation of this representation based on a period-energy 
argument [11] for Hamiltonian systems, and some properties of the periodic averaging on 
manifolds O El [T8] . This formalism allows us to get an efficient symbolic implementation 
for some models related to polynomial perturbations of the harmonic oscillator with 1:1 
resonance. In particular, we compute the second-order normal form of the Henon-Heiles [6], 
and the spring pendulum [UGH [10 J Hamiltonians, expressed in terms of the Hopf variables. 

Let us remark that the second-order normal form plays a very important role in the 
approximation of a perturbed dynamics by solutions of the averaged system when a long- 
time scale is used [21 [19]. Our desire to study this kind of dynamics led to the present 
work. 

Also, we present a package, called pdynamics, written in the CAS Maxima, which can 
automatically compute the second-order normal form in most cases of interest (see comments 
in Part [IT] for an overview of its limitations). We have chosen this particular CAS because 
of its ease of use, its syntax (very similar to that used on a blackboard), and its open-source 
character. The second part of this paper, in which we show how to use the package, is more 
elementary in mathematical terms. We give a complete list of the functions contained in 
the package with examples of use for each one of them. The software can be downloaded 
from http : / / galia . f c . uaslp . mx/~j valle j o/pdynamics . zip 
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Part I 
Theoretical background 



2 Vector fields with periodic flow 

Throughout the paper, we set S 1 = 1R/2-7TZ. We collect here some results regarding the 
flow Fl x of a vector field X, on an arbitrary manifold M, in the case when Fl x is periodic. 
Although these results are general, later they will be applied to the case of a Hamiltonian 
vector field on a Poisson manifold (M, P) . 

Let X € X{M) be a complete vector field whose flow is periodic with period function 
T € C°°(M), T > 0, that is: for any p G M, 

Fl t + T (p) = Fl x (p). (1) 



Then, X determines an S 1 — action S 1 x M — > M given by (t,p) h-> Fl^ w ^(p), where 
cj := 2ir/T > is the frequency function, and t G § . Thus, the S 1 — action is periodic, with 
constant period 2ir. 

The generator T of this S 1 — action can be readily computed: 



x v ' u)(p) ds 



t=o 

so T = -jjf. Notice, from ([TJ, that T(p) > is the period of the integral curve of X passing 
through p € M at t = 0, c p : R — >• M (which is such that c(0) = p and c p (0) = X(p)). In 
other words, c p (0) = p = c p (T(p)). Also, each point on the image of the integral curve c p , 
gives the same value for the period: T{p) = T(c p (t)), for all tel. In terms of the flow of 
X, that means 

((Fl* x )*T)(p) = T(Fl* x (p)) = T(p), for all p E M. 
As T is constant along the orbits of X, its Lie derivative with respect to X vanishes: 

d 



C * T at 



(f4)*t = 0. 

t=o 

Now, from Toj = 2ir, we get 

= C x (Tuj) = {C x T)uj + TC x oj = TC x u. 

But T > 0, so this implies that w is a first integral (or invariant) of X, 

C x u = 0. (2) 

Definition 1. A smooth function f G C°°(M) is said to be an S 1 — invariant if it is invariant 
under the flow of the generator T = ^X , that is, 

c T f = o. 

Clearly, this is equivalent to the condition (Fly)*/ = /, for all t G [0, 2n]. 
Remark 1. By ([2]), the frequency function is also an invariant of the S 1 — action: 

Cyoj = —£x^ = 0. 
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3 Averaging operators 

Given a vector field X G X(M) with periodic flow, the associated S 1 — action can be used 
to define two averaging operators, which we will denote by (•) and S. In this section, M 
will be an arbitrary manifold. 

Definition 2. For any tensor field R G VT^(M) (r—covariant, s—contravariant), the aver- 
age of R with respect to the S 1 — action on M induced by X, is the tensor field (of the same 
type as R) defined by 



(R) :-- 



2tt 



2- 



(Fr T )*M. 



The properties of the flow pQ guarantee that (R) is well-defined as a differentiable tensor 
field. Also, note that if R G YT°(M), and Xi,...,X r G X(M), a u ...,a s G ^{M) 
are arbitrary, then, for every p G M, t i-> (Fl^)* R(Xi, ... ,X r ,ai, ... ,a s )(p) is a real 
differentiable funcion on the compact [0, 2n], hence integrable. We will use this definition 
mainly applied to the case of functions / G C°°(M) ((0, 0)— tensors) and vector fields Y G 
X{M) ((0, l)-tensors). 

The other averaging operator that will be important in what follows, is the S operator. 
Definition 3. The operator S : VT s r (M) ->■ TT°(M) is given by 

S(R) :=^^ 27r (i-7r)(Fl* x )*i?dt. 

Note that both, (•) and S, are M— linear operators. Other properties are listed below. 

Lemma 3.1. For any complete vector field Y G X{M) (whose flow is not necessarily 
periodic) and smooth tensor field R G VT^(M), we have: 



_d 



s=0 



my(R) = ±-((Fi?yR-R), 



where the averaging is taken with respect to the flow ofY, that is, (R) is given by (R) := 

^CiFi^rRdt. 

Proof. Start from the identities (which follow directly from the definitions of flow and Lie 
derivative) : 



(Fl t Y y(£ Y R) = ±(Fl t Y yR= ± 



(Fl 



Y 



R 



s=0 



d_ 

d^ 



(Fl s Y )*(Fl Y )*R. 



s=0 



Taking the integral with respect to t between and 2ir on both sides, we get, on the one 
hand: 



1 



2- 



- (Fl Y nCyR)dt=- 



o 



2tt 

and, on the other: 
1 

2^ 



ds 



s=0 



1 



2tt 



2- 



(Fly)* - (FlyYRdt =- 



d 



ds 



(FlyT(R), 



2k 



(Fly)*{£ Y R) dt 



o 



2tt 



2tt 



-(Fl Y )*Rdt = — ((Ft?)* R-R). 



□ 
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Proposition 3.2. For every R G VT^(M), the following properties hold: 

(a) R is invariant under the flow of T (that is, S 1 — invariant) if and only if (R) = R. 

(b) C r {R) = 0. 

(c) If g G C°°{M) is S 1 — invariant, then (gR) = g(R). 

(d) The averaging operator commutes with tensor contractions whenever one of the tensors 
is S 1 — invariant, that is, if S G VT^(M) is S 1 — invariant and C\ is any contraction, 
then(C l k {R®S)) = C l k {{R)®S). 

Proof, (a) If R is invariant under the flow of T, then (Fl^-)*i? = R, for all t G [0, 2ir], and 
from this it is immediate that {R) = R. Reciprocally, if (R) = R we may apply the 
preceding lemma to obtain: 



d_ 



(Fl s r )*R= — ((Fl^)* R-R) 

s=0 27T 



and from the fact that the flow of T is 2tt— periodic, 

(Fl* T )*i? = 0. 



(b) From the properties of the Lie derivative and the definition of (R): 
Now, because Fl^ is 2tt— periodic: 

2tt 



(Fl* r )*(C Y (R)) = -—J (Fl u r TRdu = 0, 



so, as Fl^ is a diffeomorphism, £r(R) = 0. 

(c) It is a straightforward computation. 

(d) It is just a consequence of the commutativity between the pull-back and the tensor 
contractions, and the functorial property (Fl^)*^ <g) S) = (Fl^yR <g) (Fi^yS. 

□ 

Remark 2. In particular, from (jd]) we get that ifY G X(M) and a G J7(M) is S 1 —invariant, 
then (iyct) = i(y\a. 

Proposition 3.3. For any R G VT^(M) and g G C°°(M) S 1 — invariant, the following hold: 

(a) S(gR)=gS(R). 

(b) (C r oS)(R) = R-(R). 

Proof, (a) A straightforward computation. 
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(b) With an obvious change of variable, we have: 

{F\ s r )*S{R) = — f 27T (t - 7r)(Fl£+*)*i?dt = — (u - s - ir)(FV$)*Rdu. 

27r J 2ir J s 

Differentiating both sides of this identity with respect to the parameter s, and taking 
into account the 2ir— periodicity of the flow Fly, it results: 

The statement follows by recalling that Fl^~ is a diffeomorphism, and the identity (see 
HI): 

±(FIWS(R) = (Fl s r )*(£?S(R)). 

□ 

Finally, let us give some useful properties involving the averaging operators. 
Proposition 3.4. For all R G VT^(M), the operators Cy, {■), and S, satisfy the relations: 

(a) (C r R) = C r (R) = 0. 

(b) (S(R))=S({R}) = 0. 

(c) (da) = d(a), for all a £ U(M). 

Proof. Straightforward computations, making use of Proposition 13.31 and the fact that d 
commutes with pull-backs. □ 

4 The Hamiltonian case 

Let (M, P) be an m— dimensional Poisson manifold, where P € VA 2 TM is a Poisson bivector 
determining a bracket {/, g} = P(df, dg), for all f,g€ C°°(M). For every /, its Hamiltonian 
vector field Xf G X(M) is given by Xf(g) := {f,g}, for any g G C°°(M), equivalently, 

X f = i df P. (3) 

At any point the distribution spanned by the Hamiltonian vector fields is involutive, as 
a consequence of Jacobi's identity for the Poisson bracket {•, •}. Thus, these Hamiltonian 
vector fields give rise to a foliation whose leaves turn out to be symplectic manifolds (see 
[20]). On each leaf S, the restriction P\g is a non-degenerate Poisson bivector field which 
determines a symplectic structure as through: 

a s (X f ,X g ) :={f,g}. 

Indeed, by a theorem of Weinstein (|20j). the local structure of (M, P) can be described as 
follows: for any p G M there exists a chart (U, <f>) of M around p such that, if the coordinates 
of (/) : U -> W n (2k + I = m) are {q x , ...,q k ,pi, ...,p k ,yi, -,yi}, then 

„, d d 1 n d d 

'•^Sis'is^?;"'" < 4 » 

% — 1 — 1 



6 



where ip : ixi{U) C R l -> R is smooth and <py(p) = (717 : R m = M 2fe xlUl 1 is the 
canonical projection). The non- negative integer k is called the rank of the Poisson structure 
P at p E M. When k = m, P induces a symplectic structure on M. Then, the symplectic 
leaf S through p E M, is given by the equations (yi, ■■■,yi) = (0, ...,0). 

When moving along the flow of a Hamiltonian vector field, which is tangent to some 
integral submanifold S, it is clear that we stay on the same symplectic leaf S. Next, we 
study what happens on these leaves when the Hamiltonian vector field has periodic flow. 

We will need first an auxiliary result, interesting in its own, known as the period-energy 
relation (see [TT]). 

Proposition 4.1. Let X be a vector field on the symplectic manifold (S,o~) whose flow is 
periodic with period function T E C°°{M), T > (and frequency oj = 2tt/T). If X is the 
Hamiltonian vector field for a certain function f E C°°(M) (that is, ixo~ = —df), then: 

du A df = = dT A df. (5) 

Proof. By hypothesis, we have: 

Cx<? = ixda + dixo~ = —d 2 f = 0. 
On the other hand, using the generator T = X/oj of the S 1 — action induced by X: 

Cxo~ = uC^a do; A df. 

Recalling that to, f are first integrals of X, and hence S 1 — invariants, applying the averaging 
operator (•) to the last identity, taking into account that (Cro~) = (Proposition 13.41 (jaj)). 
and the commutativity between d and (•) (Proposition 13.41 (jcj)), we get: 

= (ujCya) - (-du A df) = oj(Cycr) - -doj Ad/ = —du A df. 

OJ ' OJ OJ 

□ 

Remark 3. Note that in terms of Hamiltonian vector fields, we can write the energy-period 
relation ([5]) as follows, 

Xaj A Xh = 0. (6) 

Also, in the course of the proof we have seen that, if T = ^X is the generator of the 
S 1 — action induced by X: 

= Cxo~ = ojCyo doj A df, 

OJ 

so from ([5]) we get the following consequence. 

Corollary 4.2. The symplectic form a is S 1 — invariant, Cya = 0. In particular, (a) = a. 

Proposition 4.3. If g E C°°(S) is 8 1 — invariant, then its Hamiltonian vector field X g E 
X(S) is also S 1 — invariant. 
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Proof. Recalling that d commutes with the averaging (Proposition 13.41 (jcj)). Remark [21 and 
the preceding Corollary, we get: 

ix ( 9 ) a = = ~^9) = (ix g <r) = i{x g )<r- 

Hence, by the non-degeneracy of a, (X g ) = X/ g \. Now, if g is 8 1 — invariant, (g) = g, and so 
(X 9 )=X g . ■ □ 

Corollary 4.4. For any S 1 — invariant g 6 C°°(S), we have 

C r X g = [X T ,X g }=0. 

Now, suppose that we are given a function H € C°°(M) on the Poisson manifold (M, P) 
such that its Hamiltonian vector field Xh £ X(M ) has periodic flow (with frequency func- 
tion uj € C°°(M), w > 0). Let T = ^Xh be the generator of the associated S 1 — action. 
From the results above we know that M is foliated by symplectic leaves S in such a way 
that P\s is equivalent to a symplectic form as (recall ([!])), and these are invariant under 
Hamiltonian flows. Thus: 

'f 1 

= C Xh P = C^P = uj£ r P A i dLU P = ujCyP + -X H A X u , 

UJ LO 

where we have used the formula C-fxA = f£xA — X A (valid for any function / € 
C°°(M), vector field X G Af(Af) and multivector field A € V(ATM), see [T?J, p. 358), as 
well as ([3]) and the fact that uj > 0. From this identity and the energy-period relation (|6|), 
we deduce that P is S 1 — invariant, 

C T P = 0. 

Moreover, if g £ C°°(M) is S 1 — invariant, the flow of its Hamiltonian vector field X g leaves 
the integral submanifolds S invariant and, as we have seen, on each of them it satisfies 
C"fX g = 0, so this is also true on M. In other words, the flows of T and X g commute on 
M. The following result exploits this fact. 

Proposition 4.5. Let (M,P) be a Poisson manifold, and H € C°°{M) such that its Hamil- 
tonian vector field X h € X(M) has periodic flow. Iff,g G C°°(M) and g is invariant under 
the induced S 1 — action, then: 

(a) {H,g}=0. 

(b) ({f,g}) = {(f),g}. 

Proof. Item (jaj) is proved by a straightforward computation. Item ([b]) is a direct consequence 
of the S 1 — invariance of g and the fact that the flows of T and X g commute. □ 

5 The main result 

Let H e = H Q +eH l + \e 2 H2+0{e i ) an e— dependent Hamiltonian function which describes a 
perturbed Hamiltonian system on a Poisson manifold (M, P), with associated bracket {•,•}. 
We will denote by Xjj e = Xjj + sXh x + \e 2 Xu 2 + 0(e 3 ) the corresponding Hamiltonian 
vector field. 
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Definition 4. The perturbed Hamiltonian vector field Xjj e is in (Deprit) normal form 
relative to Xh of order k in e if 

[X Ho ,X Hi ]=Q, for alii E {1,2,..., k}. (7) 

In terms of Hamiltonian functions, ([7]) is satisfied whenever 

{H , Hi} = 0, for all i € {1, 2, k}. 

Theorem 5.1. Suppose that the flow of Xh is periodic with frequency function uj € 
C°°(M), oj > 0. Then, there exists a canonical near-to-identity transformation <3? e on (M,a) 
which brings the perturbed Hamiltonian H £ to a Hamiltonian normal form relative to H$ of 
arbitrary order in e. In particular, the second order normal form can be expressed as: 



H £ o<S> £ = H + e(Hi) + j ((H 2 ) + ({S (^f) , Fx}}) + 0(e 3 



(8) 



Proof. If the Hamiltonian vector field Xh has periodic flow, the existence of the near-to- 
identity canonical transformation 4> e is a well known fact [31 El [T5l I16j . Here we give a 
explicit formula for it. 

Let $ £ be the flow of the perturbed vector field Z £ = Zq+eZ\ where Zq and Z\ are the Hamil- 
tonian vector field of the functions G = ^(#1) and d = i5(^ 2 + {5(^i),ili + (ili)}), 
respectively. Using the Lie transform method [BJ [U [T2l 113] . the second order development 
of H £ o $ e is given by: 

H £ o$ £ = H + e (C Zo H + H x ) 

+ £ - (C 2 Zo H + 2£ Zo Fx + C Zl Ho + H 2 ) + 0(e 3 ) (9) 

Now, we apply the results of the preceding sections to put this Hamiltonian in the form ([8]). 
To this end, we compute: 

£z Ho = -£x H S(—Hi) = (Hi) - Hi, 

CO 

L\Ho = £x Go ((Hi) - Hi) = {-S(Hi), (H t ) - H x }, 
£z Hi = Lx G Hi = {—S(Hi),Hi}, 

u UJ 

and, finally 

C Zl H = -C Xh S{H 2 + {S{-Hi),Hi + (H^}) 

u UJ 

= (H 2 ) + {{S(-Hi),Hi}) - (H 2 + {S(-Hi), Hi + {H^}). 
Substituting these identities into ([9]) , we obtain the second-order normal form (JSj) . □ 
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Further simplification in the expression of the normal form can be achieved through the 
use of Hopf variables, see subsection 19.41 below. 

Part II 
Software implementation 

To use the package, simply unpack it and copy the file pdynamics .mac in the contrib di- 
rectory of your Maxima installation (for instance, in a Unix box this directory will be some- 
thing as /usr/ share/maxima/5 . 28 . 0/share/contrib/, while in a Windows machine it will 
be located at C : \Program Files\Maxima-5 . 28 . 0\share\maxima\5 . 28 . O\share\contrib, 
depending on the version number). If you place it somewhere else, you can use Maxima's 
f ile_search command, as in f ile_search( "/home/johndoe/maxima/pdynamics .mac" ) . 
In what follows, we assume that any one of these methods has been applied, so the package 
is available to Maxima. 

(°/ il) load(pdynamics) $ 

In this part, we discuss the functions included in the package with examples of their 
usage. In the next section we offer an example of the code, the function pbracket, making 
some comments about its implementation, but due to space reasons we omit the code for 
the remaining functions. The interested reader can take a look at the source of the package, 
which for the most part is self-explanatory. 

Let us say some words about the limitations of the package. As it is based on the 
computation of averages along the flow of a Hamiltonian vector field, it involves at some 
point the symbolic computation of an integral (when determining the integral curves that 
define the flow). Thus, the package is as good as it is the Maxima symbolic integrator, 
which is far from perfect. It may happen that some integral can be done "by hand" and 
Maxima can not solve it, or that some other CAS (like Maple or Mathematica) can find 
the integral while Maxima can not. In these cases, the user can circumvent the difficulty by 
using this external output to directly define the Hamiltonian flow of the Hamiltonian under 
study, say h, as phamflow(h) : = [F1, . . . ,Fn] , where Fl,...,Fn are the components of the 
flow, found by whatever means. But, of course, when some particular integral not solvable 
in closed form appears, the package is useless and numerical methods are required. 

Another drawback is that the frequency function u must be supplied by the user. Even 
in the simplest cases, to determine in an automated way whether a certain function is 
periodic and, if so, to compute its (shortest) period is a tricky task (for instance, consider 
the Dirichlet function Xq), for Maxima or for any other CAS3- Thus, we have preferred 
to avoid it here. In most practical cases, the Hamiltonian will be a perturbation of the 
harmonic oscillator, so this seems to be not a serious problem^. 

x For a discussion related to Wolfram Alpha, see: 

: //math, stackexchange . com/questions/ 141033/how-to-ef feet ively- compute- a-per iodic- function 
2 Here is an example of what the user can do in most cases. Suppose we want to find the period of 
the function fit) — pi sin u)t / 'uj + q\ cosuit (see 18.71 below). In a Maxima session, do load(to_poly_solve) , 
then f (t) : =pl*sin(w*t) /w+ql*cos (w*t) , and nicedummies C/,solve(f (t+T) -f (t) ,T) ) . Among the solu- 
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6 Poisson brackets 



Here we implement a function for computing the Poisson bracket of two functions f,g& K : 

If 1 - V ( d f dg d f d 9 \ 
^ \dpi dqi dqi dpi] ' 

We must define first the functions in the form f(qi,pi,...,q n ,p n ). The order of the arguments 
is important, but their names are not. The function pbracket always returns the answer 
in the standard form (with coordinates qi,pi, ...,q n ,p n ). 

pbracket (f , g) : =block ( [n , Q , P , vars] , 

n : length (args(lhs (apply (fundef, [f] ))))/2, 

Q rmakelist (concat (q,j),j,l,n), 

P :makelist (concat (p, j) , j , l,n) , 

vars: join(Q,P) , 

lsum( 

-dif f (apply (f , vars) , Q [i] ) *dif f (apply (g, vars) ,P [i] ) 
+dif f (apply (f , vars) ,P [i] ) *dif f (apply (g, vars) , Q [i] ) , 
i ,makelist ( j , j , 1 ,n) ) 
)$ 



The first lines of this function appear repeatedly in other functions, so let us briefly 
explain what they do. As the dimension n in M. 2n is variable, the first thing to do is to 
know it, and this can be achieved by looking at the number of arguments of the functions 
/, g. The second line does thatd. Thus, if / is of the form f = f(a, b, c, d) we know that 
n = 2. Next, we form a list of the variables involved, naming them internally in a consistent 
way as (ql,pl, qn,pn). To work iteratively we need a labelling and we have chosen it 
to be (qj,Pj) (note that given two lists Q : [ql, qn], P : [pi, ...,pn], Maxima's command 
join will intersperse them). This labelling is done independently of the user's one, so if the 
functions have been defined as, say, / = f(a, b, c, d) and g = g(x, y, z, t), we will treat them 
as functions / = f(ql,pl,q2,p2),g = g(ql,pl,q2,p2) and the output {/,<?} will be again a 
function of (ql,pl, q2,p2). The user can change then the name of the variables to whatever 
she wants. 

6.1 Example 

Consider in M 4 with coordinates (qi,Pi, <72,P2) (we denote q = ((71,(72)1 P = (pi,P2) and 
\q\ = V<li + 121 bl = yPl +P2) tne functions ^\q\ 2 , ^\p\ 2 and q ■ p = q\p\ + q 2 P2 ■ 

(%i2) normq(ql,pl,q2,p2) : =(ql~2+q2~2) /2$ 
(%i3) normp(ql,pl,q2,p2) : =(pl~2+p2"2) /2$ 



tions, there is the one of interest for us: 2irk/uj, k £ Z. 

3 Thanks to R. Dodier, from the Maxima list, for this trick. 
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C/„i4) prodqp(ql,pl,q2,p2) : =ql*pl+q2*p2$ 



These functions close in (C°°(IR 4 ), {.,.}) a Lie subalgebra isomorphic to sla(K): 
(°/ i5) pbracket(normp,normq) ; 

(%o5) p2 q2+plql 

(°/ i6) pbracket(prodqp,normq) ; 

(%o6) q2 2 + ql 2 

(°/ i7) pbracket(prodqp,normp) ; 
(%o7) - P 2 2 - pi 2 
6.2 Example 

Consider a central force in M 3 , i 7 , with potential V = U(r), where r = \J x 2 + y 2 + z 2 . The 
components of F are those of the gradient — gradV, and they are computed using the chain 
rule, as the dependence of U with r is undetermined. We can work in Maxima with an 
arbitrary U(r) declaring its (also undetermined) derivative to be U'(r): 

(°/ i8) V(x,y,z) :=U(sqrt(x~2+y~2+z~2))$ 
C/.i9) gradef(U(r),U\'(r))$ 

For instance, F x = —dV/dx is given by 
C/.ilO) diff (V(x,y,z) ,x) ; 



(%o!0) 




\J z 2 + y 2 + x 2 

For a particle moving in M 3 under the influence of U, its Hamiltonian is 



(°/ ill) Hcentral(x,px,y,py,z,pz) : = (px~2+py"2+pz~2)/2 + V(x,y,z)$ 



and its angular momentum in the z direction 



(%il2) L[z] (x.px.y.py.z.pz) :=x*py-y*px$ 
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Then, no matter what the explicit expression of U is, the z component of the momentum 
is an integral of motion: 

(°/ il3) pbracket(Hcentral,L[z] ) ; 

(%ol3) 

The same is true for the other components of L = (L x , L y , L z f?\. 

7 Hamilton's equations 

Given a Hamiltonian H(qi,pi, ...,q n ,p n ), the evolution equation for any classical observable 
/ = f(Ql,Pl,—,Qn,Pn) along a physical trajectory c(i) = (qi(t),pi(t), q n (t),p n (t)) in 
phase space, is 

^ = {",/}• (10) 

In particular, the canonical (Hamilton's) evolution equations are: 

dq { _ dH 
dt dpi 

dt % ' 1 j 

for i £ {1, ...,n}. In order to write down the evolution equation for an observable, we first 
solve (fTTj) and then substitute the obtained solutions qi = qi(t),pi = Pi(t) into ((TOj) . 

7.1 Canonical equations 

The following functions return and solve the canonical equations. The Hamiltonian must 
be defined first in the form H = H{q\,p\, q n ,p n )- The ordering of the arguments of H 
is important, but their names are not. The function pcanonical_eqs gives a list of the 
form [eq gi , eq Pl , eq qn , eq Pn ] where each eq qi , eq Pi is the first-order equation corresponding 
to qi(t), Pi(t), respectively. On the other hand, pcanonical_sol returns a list of the form 
[qi(t),pi(t), q n (t),p n (t)], the solutions to these equations. 

7.2 Example 

The freely falling particle in an homogeneous gravitational field has a Hamiltonian: 
(%il4) K(q,p) :=p~2/2+m*g*q$ 



so the canonical equations are: 

4 In fact, it is true that the square norm of the total angular momentum L 2 = l? x + Ly + L 2 Z has vanishing 
Poisson bracket with H central, but we can't compute directly pbracket (Hcentral , L 2 . + Ly + L 2 ) because 
L 2 +L 2 + L 2 is not a single explicit function. While it is possible to implement properties such as bilinearity, 
the Jacobi identity or the Leibniz rule in the definition of pbracket, the form presented here is sufficient for 
our purposes. 
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(°/ il5) pcanonical_eqs (K) ; 



(%ol5) ql (t) = pi (t) , A pi (t) = -g m] 
with solutions: 

(°/ il6) pcanonical_sol(K) ; 

Q Tfl t 

(%ol6) [ql (t) = - + pi (0) t + ql (0) , pi (t) = pi (0) - m t] 

7.3 Example 

For the harmonic oscillator of frequency uj > 0, with Hamiltonian: 
(°/ il7) assume (%omega > 0) ; 

(%ol7) [u > 0] 

(%il8) Hosc(x,y) :=y~2/2+°/.omega~2*x~2/2; 

(%ol8) Rosc(x,y) ;=t- + ^- 
the canonical equations are: 
(°/ il9) pcanonical_eqs (Hose) ; 

(%ol9) [j- ql (t) = pi (t) , pi (t) = ql (t)] 
at at 

with solutions: 

(°/„i20) pcanonical_sol(Hosc) ; 

(%o20) [ql(t) = P 1 ^) sm M) +ql (Q) C os(wt) ,pl (t) =pl(0) cos (wt)-ql (0) wsin(wt)] 

7.4 Evolution of observables 

The function pevolution computes the evolution of an observable / = f(qi,Pi,---,q n ,Pn) 
with respect to a Hamiltonian H = H(qi,pi, q n ,Pn)- It needs that the functions pbracket 
and pcanonical_sol be previously loaded. 

7.5 Example 

Let us do a little sanity check. If the observable is one of the coordinates qi, its evolution 
must coincide with that resulting from the canonical equations. In this case we compute 
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the evolution of the ql coordinate of the harmonic oscillator (so the result should coincide 
with the first component in (%o20)): 

(%i21) foo(a,b) :=a$ 

(°/ i22) pevolution(f oo,Hosc) ; 

(%o22) F (t) = Pl(0) S ' m{UJt) + ql (0) cos (u, t) 

8 Vector fields and flows 

8.1 Integral curves of a vector field 

The function pintcurv computes the integral curves of a vector field X in R m (not neces- 
sarily even dimensional) given by the list of its components [Xi, ...,X m ], each one of them 
a function of the local coordinates [xi, x m }. The user must supply the list of components 
of X and a list containing the names of the coordinates. 

8.2 Example 

The integral curves of the vector field X = (ax — b, ay) in the plane can be computed as 
follows: 

(°/ i23) pintcurv ( [a*x-b,a*y] , [x,y]); 

(%o23 ) [x(t ) = i - "-»•"" ,,(!) = y(0) e -] 
a a 

8.3 The flow of a vector field 

If X is a vector field on M m , its flow is a mapping Fix : K x M. m — s- M. m such that, if c p (t) is 
the integral curve passing by p at t = 0, 

Fl x (t,p) = F\ t x (p) = c p (t). 

The function pvectflow computes the flow of the vector field X = (X\, X m ). As in 
the previous case, the input is a couple of lists: one containing the components of X, 
[X\, ...,X m ], and the other containing the coordinates used on M m , [x%, ...,x m ]. It returns 
a list with the flow mapping components [(Flx)i(t,xi, ...,x m ), ...,(Flx) m (t,xi, ...,x m )]. 

8.4 Example 

The flow of the vector field X = (y, —lo 2 x) in the plane is: 
(°/ i24) pvectf low( [y,-°/ omega~2*x] , [x,y]); 

sin (cot) y 

(%o24) [ h cos [lo t) x, cos (uit) y — u sm [u t) x\ 
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8.5 The Hamiltonian vector field 



Assume the manifold IR endowed with the canonical symplectic form 



f2 = dpi A dqi + • • • + dp n A dq. 



(12) 



Given a Hamiltonian H = H (qi,p±, q n ,Pn), its associated Hamiltonian vector field Xjj 
is given by the condition ix H Q = —dH. This is easily seen to lead to the components 



in the basis {d/dqi, d/dpi, d/dq n , d/dp n }. The function phamvect computes Xh from 
H, expressing its components in the form [X\(q,p), X2 n (q,p)]- 

8.6 Example 

For the harmonic oscillator, the Hamiltonian vector field is: 
(%i25) phamvect (Hose) ; 

(%o25) \pl,-uj 2 ql] 

8.7 The Hamiltonian flow 

Suppose we have a Hamiltonian vector field on R 2n endowed with the canonical symplectic 
form (fT2|) above. Given a Hamiltonian H = H(qi,p\, q n ,p n ), the flow of its associated 
vector field Xh is called the Hamiltonian flow. 

The function phamf low computes the Hamiltonian flow determined by a Hamiltonian H. It 
is just the composition of pvectf low and phamvect. 

For example, if H is taken to be the Hamiltonian of the harmonic oscillator we recover 
previous results (cfr. Examples 18.41 and I8.6P : 

(°/ i26) phamf low (Hose) ; 



9 The averaging method for normal forms 
9.1 Averaging of a function respect to a periodic flow 

Suppose we have a Hamiltonian system (with phase space W 11 ) on which there is an 
S 1 — action with generator X. Then, the flow Fl^ is periodic. The average of an observable 
g with respect to the induced S 1 — action is the function defined by 




(%o26) [ 



pi sin (oj t) 



+ ql cos (a; t) , pi cos (wt) -w ql sin (w t)] 
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The paverage function below computes (g) given the input g, X = [Xi, ...,X m ] (a 
list containing the components of the vector field X in the basis {d/dxi, ...,d/dx m }) and 
x = [xi, ...,x m ] (a list with the coordinates used in M. m ). The function g must have been 
previously defined. 

Usually, the action is Hamiltonian, that is, the vector field X = Xh for some Hamilto- 
nian H. Assuming in this case that the phase space is M. 2n with the canonical symplectic 
form (|12p . the function phamaverage computes the average of g with respect to the Hamil- 
tonian vector field of H. 

(°/ i27) goo(x,y) : =-x~2* (1+y) /2$ 
(°/ i28) paverage (goo , [y , -x] , [x , y] ) ; 

(%o28) - V —^- 

(°/ i29) HoscO(ql,pl) : =(ql~2+pl~2) /2$ 
(%i30) phamaverage (goo, HoscO) ; 

(%o30) 



ql 2 + pl 2 



4 

There is another average which is very important in the theory of normal forms. It is 
given by the action of the operator S: 



S(g) = ±J*\t-n)(Fl t x ygdt. 



The command psprojector computes it: 
(°/ i31) psprojector (goo, [y,-x] , [x,y]); 

(% o3 l) 3x ^~ 2x3 

9.2 Second-order normal form of a perturbed Hamiltonian 

The previous routines are all we need for computing the normal form of a Hamiltonian on 
M 2 endowed with the canonical symplectic form (|12p . If we have a system admitting an 
S 1 — action, described by a perturbed Hamiltonian 

H = H + eH 1 + ^H 2 , 

and such that the Hamiltonian vector field of Hq, Xh , has periodic flow with frequency u), 
then its second-order normal form is given by 
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The command pnormal2 performs the necessary computations given the Hamiltonians Hq, 
Hi, H2, and the parameter e. Another variable uj (the frequency function for the flow of 
Xh ) is optional: if it is not included, it is assumed that this frequency is uj = 1. That 
function uj, if included in the argument of pnormal2, must have been previously defined. 

9.3 Example: The Henon-Heiles Hamiltonian 

This example is taken from [7j. The Hamiltonian is 

K= 1 -(pl+p 2 2 )+ 1 -( q 2 l +q 2 2 )+e( q 4- qiq j 



(note that the perturbation term is an homogeneous polynomial of degree 3), so we define: 
(°/ i32) K0(ql,pl,q2,p2) : = (pl~2+p2~2)/2+(q:r2+q2~2)/2; 

(%o32) K0(ql,pl,q2,p2) := P - ^ + ^±^! 
(%i33) Kl(ql,pl,q2,p2) : =ql~3/3-ql*q2~2 ; 

(%o33) Kl(ql,pl,q2,p2) : =^-ql q 2 2 
(7.134) K2(ql ) pl ) q2,p2) :=0; 

(%o34) K2(ql,pl,q2,p2) := 

The frequency function for the flow of Xk is readily found to be (see footnote [2] in page 

ED: 

C/.135) u(ql,pl,q2,p2):=l$ 

The second-order normal form is therU: 
(%i36) pnormal2(K0 ) Kl,K2 ) %epsilon,u) ; 

(%o36) 

+ " ^ (5^ + (10 gl 2 + 10p2 2 - 18pl 2 ) ,2 2 

+56plp2gl q2 + 5ql 4 + (lOpl 2 - 18p2 2 ) ql 2 + 5p2 4 + Wpl 2 p2 2 + 5pl 4 ) 

9.4 Hopf variables 

It is usual to express the normal form in terms of the Hopf variables w\, w%, W3, w±, as a 
previous step to carry on the reduction of symmetry process (see [6],[7j). For the case in 
which Hq is the Hamiltonian of the 2D— harmonic oscillator, these variables form a system 



We have slightly edited the output in order to make it more readable. 
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of functionally independent generators of the algebra of first integrals of Hq, and are denned 
as w\ = 2(q 1 q 2 +pip 2 ), w 2 = 2(q 1 p 2 -q2Pi), w 3 = ql+p\-ql~p% = ql + ql+p\+pl- The 
functions phopf 2, phopf 4 attempt to express a given expression (a homogeneous polynomial 
in the variables (qi, q 2 ,pi,p 2 ) of degree 2 and 4, respectively) in terms of them. To apply 
these functions to the output of pnormal2 above, we can select the independent term and 
the coefficient of e 2 as follows: 

(%i37) phopf2(coeff (°/ 0) °/ epsilon,0) ) ; 

(%c37) S 

(%i38) phopf 4(coeff (7„th(2) ,%epsilon"2) ) ; 

,„ OQ , w 2 2 (48%rl + 7) ^i(48%rl + 5) 2 2 

(%o38) — ^ — — — + wi %rl + wi %rl 

48 48 6 

The formulas appearing in [7] are recovered by choosing the value of the parameter: 
C/.i39) subst(%rl=0,%) ; 

Thus, the second-order normal form of the Henon-Heiles system is 
H t o^ = ^ + ^(7 W l-5 W i)+0(e 3 ). 
9.5 Example: The spring pendulum 

Consider the case of the Hamiltonian of a spring-pendulum (sec [5j,[4 ,[10j): 

H(ql,pl, q2,p2) = + - \q\{l + q 2 ), 

which is that of a perturbed system Hq + eH\, where 

Ho(5 i, P i, ? 2, P 2) = ^ + ii±i, 

and 

Hi(ql,pl,q2,p2) = . 

Note that the perturbation term now is not homogeneous. We define the terms of the 
perturbed Hamiltonian: 

(%i40) H0(ql,pl,q2,p2) :=(pl~2+p2~2)/2+(ql-2+q2~2)/2; 

(%o40) HO (ql,pl,q2,p2) := V — % V — + Q — ± ^ 
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C/.i41) Hl(ql,pl,q2,p2) :=-ql~2*(l+q2)/2; 



(%o41) m(ql,pl, q 2,p2) := ( ^ ^ + g2) 
C/.142) H2(ql,pl,q2,p2) :=0; 

(%o42) H2(gl,pl,g2,p2) := 

and compute the normal form in the orig inal variable^. Note that we do not explicitly 
write the frequency function (thus assuming it is the constant 1): 

C/.i43) pnormal2 (HO , HI , H2 , %epsilon) ; 

(%o43) 

P 2* +P l* q2* + q l*_e 

2 2 4 v ; 

~ 192 (( 20ql2 ~ A P l2 ) <?2 2 + 48plp2glg2 + 5gl 4 + 
(-4p2 2 + 10pl 2 + 12) ql 2 + 20pl 2 p2 2 + 5pl 4 + 12pl 2 ) 

As before, we can express in terms of the Hopf variables the independent terms and the 
coefficient of e: 

(%i44) phopf2(coeff (°/ 0) °/ epsilon,0) ) ; 

W4 

~2 

(%i45) phopf2(coeff (7„th(2) ,%epsilon, 1) ) ; 



(%o44) 2 



(%o45) -f-f 

Note that the coefficient of e 2 is not a homogeneous polynomial (of degree 4): there are 
two 2— degree terms: (qf + p 2 )/16. Thus, it does not make sense to apply phopf4, as this 
would lead to an error. Luckily, these terms can be easily expressed in terms of the variables 
w%, W2-, Ws, W4 (as {q\ + p 2 )/16 = (u>4 + wz)/32) and then we can analyse the remainder, 
which is a polynomial of degree 4: 

(%i46) phopf4(coeff (%th(3) ,%epsilon,2)+12*(ql~2+pl~2)/192) ; 

tv a#\ ^!(768%r2 + 25) w 2 3 (256 %r2 + 5) «; 2 (32%r2 + l) 2 5u^ 4 

( /oo4o I h tf 1 /orz 

v ' 768 256 32 1 384 

Let us take the simplest solution: 



Again, we have slightly edited the output. 
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C/.i47) subst(°/,r2=0,%) ; 



(%o47) - — - + — + — 

The leftover in the coefficient of e 2 is 
(%i48) phopf2(12*(ql-2+pl"2)/192) ; 



(%o48) ^ + ^ 
v ; 32 32 



Thus, we get the second-order normal form of the spring pendulum in the Hopf variables: 

wa e , e 2 / 9 25wa 2 5w^wa 5w^ 2 \ o, 
if e o $ e = -± - - (w 4 + w 3 ) + — [wa + w 3 + w\ -= -2-* + — ^ + 0(e 3 ). 



2 8 v ' 32 V 24 12 
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